Construction d’une typologie de piétons à partir d’une enquête ménages déplacements et représentation cartographique des profils : l’exemple de Bogotá

1 Calcul des profils de piétons au moyen d’une AFDM et d’une CAH

1.1 Chargement des packages nécessaires

#Chargement des packages
library(FactoMineR) # Pour realiser une analyse multivariee
library(openxlsx) # Pour importer/exporter depuis/vers Excel
library(factoextra) # pour produire des graphiques esthetiques a partir des resultats de l'AFDM
library(ggplot2) # package necessaire au fonctionnement de factoextra
library(sf) # pour manipuler les objets ayant une composante spatiale
library(spatstat) # pour le lissage spatial
library(sp) # pour utiliser la fonction as.owin de spatstat
library(maptools) # pour utiliser la fonction as.owin de spatstat
library(cartography) # pour la representation cartographique et l'habillage des cartes
library(raster) # pour le traitement de donnees matricielles
library(tidyverse) # pour utiliser la fonction "separate"

1.2 Import du jeu de données et résumé statistique

VariablesTypoUltraMarcheurs <- read.xlsx("02_analyse_ultra_marcheurs_v6.xlsx")
rownames(VariablesTypoUltraMarcheurs) <- VariablesTypoUltraMarcheurs$id_hog_pers # pour ajouter l'id des personnes comme entête de ligne
str(VariablesTypoUltraMarcheurs) # Pour connaître le type des variables et avoir un aperçu des valeurs prises pour chaque variable

## 'data.frame':    15482 obs. of  24 variables:
##  $ id_hog                       : chr  "P01283" "8856" "5268" "3478" ...
##  $ id_hog_pers                  : chr  "P01283-4" "8856-1" "5268-2" "3478-1" ...
##  $ ocasional                    : chr  "regular" "regular" "regular" "regular" ...
##  $ fq_sem_viaje                 : num  39765 30081 23469 22180 19883 ...
##  $ duracion_prom                : num  24.5 36.5 8 23.8 30 ...
##  $ calif_exp                    : num  5 4.81 3.14 5 5 ...
##  $ volv_casa                    : num  19883 15598 10058 11090 9941 ...
##  $ trabajo                      : num  0 0 6705 0 0 ...
##  $ escuela                      : num  0 0 0 0 9941 ...
##  $ cuidado                      : num  19883 10027 3353 11090 0 ...
##  $ compras                      : num  0 4456 3353 0 0 ...
##  $ edad                         : chr  "15-39" "15-39" "15-39" "15-39" ...
##  $ niv_edu                      : chr  "prim_sec" "univ" "media" "media" ...
##  $ genero                       : chr  "muj" "muj" "muj" "muj" ...
##  $ difisi                       : chr  "difisi_no" "difisi_no" "difisi_no" "difisi_no" ...
##  $ licond                       : chr  "licond_no" "licond_no" "licond_no" "licond_no" ...
##  $ estratoCL                    : chr  "s_1-2" "s_1-2" "s_1-2" "s_1-2" ...
##  $ tot_pers                     : num  7 5 4 3 7 5 3 2 5 2 ...
##  $ num_moto                     : num  0 0 1 0 0 1 0 0 0 0 ...
##  $ num_carro                    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num_bici                     : num  0 1 1 2 0 0 1 1 0 0 ...
##  $ GrandesZones                 : chr  "Sur" "Sur" "Sur" "Sur" ...
##  $ f_exp                        : num  1988 1114 671 584 1988 ...
##  $ ZAT_Dest_1er_viaje_desde_casa: num  544 714 351 993 567 663 114 115 547 452 ...

head(VariablesTypoUltraMarcheurs)

##          id_hog id_hog_pers ocasional fq_sem_viaje duracion_prom calif_exp
## P01283-4 P01283    P01283-4   regular     39765.36      24.50000  5.000000
## 8856-1     8856      8856-1   regular     30081.17      36.48148  4.814815
## 5268-2     5268      5268-2   regular     23468.56       8.00000  3.142857
## 3478-1     3478      3478-1   regular     22180.43      23.81579  5.000000
## P01283-6 P01283    P01283-6   regular     19882.68      30.00000  5.000000
## 3690-2     3690      3690-2   regular     19044.74       9.00000  5.000000
##          volv_casa  trabajo  escuela   cuidado  compras  edad  niv_edu genero
## P01283-4 19882.678    0.000    0.000 19882.678    0.000 15-39 prim_sec    muj
## 8856-1   15597.642    0.000    0.000 10027.056 4456.469 15-39     univ    muj
## 5268-2   10057.953 6705.302    0.000  3352.651 3352.651 15-39    media    muj
## 3478-1   11090.215    0.000    0.000 11090.215    0.000 15-39    media    muj
## P01283-6  9941.339    0.000 9941.339     0.000    0.000  6-15 prim_sec    muj
## 3690-2    7617.896    0.000    0.000  7617.896 3808.948 40-64 prim_sec    muj
##             difisi    licond estratoCL tot_pers num_moto num_carro num_bici
## P01283-4 difisi_no licond_no     s_1-2        7        0         0        0
## 8856-1   difisi_no licond_no     s_1-2        5        0         0        1
## 5268-2   difisi_no licond_no     s_1-2        4        1         0        1
## 3478-1   difisi_no licond_no     s_1-2        3        0         0        2
## P01283-6 difisi_no licond_no     s_1-2        7        0         0        0
## 3690-2   difisi_no licond_no     s_1-2        5        1         0        0
##          GrandesZones     f_exp ZAT_Dest_1er_viaje_desde_casa
## P01283-4          Sur 1988.2678                           544
## 8856-1            Sur 1114.1173                           714
## 5268-2            Sur  670.5302                           351
## 3478-1            Sur  583.6956                           993
## P01283-6          Sur 1988.2678                           567
## 3690-2            Sur  761.7896                           663

summary(VariablesTypoUltraMarcheurs) # Résumé des valeurs contenues dans l'objet

##     id_hog          id_hog_pers         ocasional          fq_sem_viaje    
##  Length:15482       Length:15482       Length:15482       Min.   :    0.0  
##  Class :character   Class :character   Class :character   1st Qu.:    0.0  
##  Mode  :character   Mode  :character   Mode  :character   Median :  567.6  
##                                                           Mean   : 1069.7  
##                                                           3rd Qu.: 1521.6  
##                                                           Max.   :39765.4  
##                                                                            
##  duracion_prom       calif_exp       volv_casa          trabajo       
##  Min.   :   0.00   Min.   :0.000   Min.   :    0.0   Min.   :   0.00  
##  1st Qu.:   9.00   1st Qu.:3.667   1st Qu.:   93.9   1st Qu.:   0.00  
##  Median :  15.00   Median :4.857   Median :  295.4   Median :   0.00  
##  Mean   :  21.09   Mean   :4.085   Mean   :  553.1   Mean   :  86.79  
##  3rd Qu.:  25.00   3rd Qu.:5.000   3rd Qu.:  756.5   3rd Qu.:   0.00  
##  Max.   :1080.00   Max.   :9.000   Max.   :19882.7   Max.   :7105.89  
##                                                                       
##     escuela          cuidado           compras            edad          
##  Min.   :   0.0   Min.   :    0.0   Min.   :   0.00   Length:15482      
##  1st Qu.:   0.0   1st Qu.:    0.0   1st Qu.:   0.00   Class :character  
##  Median :   0.0   Median :    0.0   Median :   0.00   Mode  :character  
##  Mean   : 221.9   Mean   :  140.7   Mean   :  92.85                     
##  3rd Qu.: 181.3   3rd Qu.:    0.0   3rd Qu.:   0.00                     
##  Max.   :9941.3   Max.   :19882.7   Max.   :5265.09                     
##                                                                         
##    niv_edu             genero             difisi             licond         
##  Length:15482       Length:15482       Length:15482       Length:15482      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##   estratoCL            tot_pers         num_moto        num_carro     
##  Length:15482       Min.   : 1.000   Min.   :0.0000   Min.   :0.0000  
##  Class :character   1st Qu.: 3.000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Mode  :character   Median : 4.000   Median :0.0000   Median :0.0000  
##                     Mean   : 3.908   Mean   :0.1723   Mean   :0.2958  
##                     3rd Qu.: 5.000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##                     Max.   :13.000   Max.   :8.0000   Max.   :6.0000  
##                                                                       
##     num_bici       GrandesZones           f_exp        
##  Min.   : 0.0000   Length:15482       Min.   :   1.00  
##  1st Qu.: 0.0000   Class :character   1st Qu.:  49.69  
##  Median : 1.0000   Mode  :character   Median : 108.10  
##  Mean   : 0.9474                      Mean   : 140.62  
##  3rd Qu.: 2.0000                      3rd Qu.: 189.25  
##  Max.   :50.0000                      Max.   :1988.27  
##                                                        
##  ZAT_Dest_1er_viaje_desde_casa
##  Min.   :   0.0               
##  1st Qu.: 351.0               
##  Median : 571.0               
##  Mean   : 615.8               
##  3rd Qu.: 794.0               
##  Max.   :1908.0               
##  NA's   :766

1.3 Apurement du jeu de données

VariablesTypoUltraMarcheurs <- na.omit(VariablesTypoUltraMarcheurs) # suppression des individus pour lesquels des valeurs ne sont pas renseignées

VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$num_bici <= 10),] # suppression des individus résidant dans des ménages dans lesquels il y a plus de 10 vélos (situation peu probable, sans doute des erreurs de saisie)

VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$calif_exp <= 5),] # suppression des individus ayant des qualifications d'experiences de trajet superieures à 5

VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$duracion_prom <= 300),] # suppression des individus ayant une durée moyenne de déplacement à pied supérieure à 5 heures (300 minutes)

1.4 Suppression des colonnes inutiles pour l’AFDM

# colonnes d'identification des individus
VariablesTypoUltraMarcheurs$id_hog <- NULL 
VariablesTypoUltraMarcheurs$id_hog_pers <- NULL

# colonnes contenant des indications de localisation spatiale
VariablesTypoUltraMarcheurs$ZAT_Dest_1er_viaje_desde_casa <- NULL 

# Autres colonnes 
VariablesTypoUltraMarcheurs$f_exp <- NULL
View(VariablesTypoUltraMarcheurs) # Resume variables
summary(VariablesTypoUltraMarcheurs)

##   ocasional          fq_sem_viaje     duracion_prom      calif_exp    
##  Length:14218       Min.   :    0.0   Min.   :  1.00   Min.   :1.000  
##  Class :character   1st Qu.:    0.0   1st Qu.: 10.00   1st Qu.:4.000  
##  Mode  :character   Median :  628.5   Median : 15.00   Median :5.000  
##                     Mean   : 1109.0   Mean   : 20.83   Mean   :4.234  
##                     3rd Qu.: 1548.6   3rd Qu.: 25.00   3rd Qu.:5.000  
##                     Max.   :39765.4   Max.   :291.25   Max.   :5.000  
##    volv_casa          trabajo           escuela          cuidado       
##  Min.   :    0.0   Min.   :   0.00   Min.   :   0.0   Min.   :    0.0  
##  1st Qu.:  111.0   1st Qu.:   0.00   1st Qu.:   0.0   1st Qu.:    0.0  
##  Median :  323.0   Median :   0.00   Median :   0.0   Median :    0.0  
##  Mean   :  574.5   Mean   :  87.98   Mean   : 233.2   Mean   :  144.2  
##  3rd Qu.:  769.5   3rd Qu.:   0.00   3rd Qu.: 223.7   3rd Qu.:    0.0  
##  Max.   :19882.7   Max.   :7105.89   Max.   :9941.3   Max.   :19882.7  
##     compras            edad             niv_edu             genero         
##  Min.   :   0.00   Length:14218       Length:14218       Length:14218      
##  1st Qu.:   0.00   Class :character   Class :character   Class :character  
##  Median :   0.00   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :  96.57                                                           
##  3rd Qu.:   3.00                                                           
##  Max.   :5265.09                                                           
##     difisi             licond           estratoCL            tot_pers     
##  Length:14218       Length:14218       Length:14218       Min.   : 1.000  
##  Class :character   Class :character   Class :character   1st Qu.: 3.000  
##  Mode  :character   Mode  :character   Mode  :character   Median : 4.000  
##                                                           Mean   : 3.886  
##                                                           3rd Qu.: 5.000  
##                                                           Max.   :13.000  
##     num_moto        num_carro        num_bici       GrandesZones      
##  Min.   :0.0000   Min.   :0.000   Min.   : 0.0000   Length:14218      
##  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.: 0.0000   Class :character  
##  Median :0.0000   Median :0.000   Median : 1.0000   Mode  :character  
##  Mean   :0.1678   Mean   :0.288   Mean   : 0.9254                     
##  3rd Qu.:0.0000   3rd Qu.:0.000   3rd Qu.: 2.0000                     
##  Max.   :8.0000   Max.   :5.000   Max.   :10.0000

str(VariablesTypoUltraMarcheurs)

## 'data.frame':    14218 obs. of  20 variables:
##  $ ocasional    : chr  "regular" "regular" "regular" "regular" ...
##  $ fq_sem_viaje : num  39765 30081 23469 22180 19883 ...
##  $ duracion_prom: num  24.5 36.5 8 23.8 30 ...
##  $ calif_exp    : num  5 4.81 3.14 5 5 ...
##  $ volv_casa    : num  19883 15598 10058 11090 9941 ...
##  $ trabajo      : num  0 0 6705 0 0 ...
##  $ escuela      : num  0 0 0 0 9941 ...
##  $ cuidado      : num  19883 10027 3353 11090 0 ...
##  $ compras      : num  0 4456 3353 0 0 ...
##  $ edad         : chr  "15-39" "15-39" "15-39" "15-39" ...
##  $ niv_edu      : chr  "prim_sec" "univ" "media" "media" ...
##  $ genero       : chr  "muj" "muj" "muj" "muj" ...
##  $ difisi       : chr  "difisi_no" "difisi_no" "difisi_no" "difisi_no" ...
##  $ licond       : chr  "licond_no" "licond_no" "licond_no" "licond_no" ...
##  $ estratoCL    : chr  "s_1-2" "s_1-2" "s_1-2" "s_1-2" ...
##  $ tot_pers     : num  7 5 4 3 7 5 3 2 5 2 ...
##  $ num_moto     : num  0 0 1 0 0 1 0 0 0 0 ...
##  $ num_carro    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num_bici     : num  0 1 1 2 0 0 1 1 0 0 ...
##  $ GrandesZones : chr  "Sur" "Sur" "Sur" "Sur" ...
##  - attr(*, "na.action")= 'omit' Named int [1:767] 77 78 185 279 288 539 594 595 596 650 ...
##   ..- attr(*, "names")= chr [1:767] "964-1" "964-3" "P00002-2" "4497-1" ...

1.5 Calcul de la première AFDM (Analyse factorielle des Données Mixtes)

res <- FAMD(VariablesTypoUltraMarcheurs, ncp = 10, graph = FALSE)

#Affichage des plans factoriels
plot(res, (choix = c("ind", "var", "quanti", "quali")), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

plot(res, choix = c("var"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes") 

plot(res, choix = c("quanti"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes") 

plot(res, choix = c("quali"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

#Affichage de la fréquence cumulée des valeurs propres
get_eigenvalue(res)

##        eigenvalue variance.percent cumulative.variance.percent
## Dim.1    3.408112        10.993911                    10.99391
## Dim.2    2.457661         7.927940                    18.92185
## Dim.3    1.871238         6.036250                    24.95810
## Dim.4    1.661313         5.359074                    30.31717
## Dim.5    1.303185         4.203824                    34.52100
## Dim.6    1.265837         4.083345                    38.60434
## Dim.7    1.120536         3.614631                    42.21897
## Dim.8    1.099534         3.546882                    45.76586
## Dim.9    1.079988         3.483832                    49.24969
## Dim.10   1.040306         3.355826                    52.60551

fviz_screeplot(res)

1.6 Affichage des variables pour en placer en supplémentaire

print(as.data.frame(colnames(VariablesTypoUltraMarcheurs)))

##    colnames(VariablesTypoUltraMarcheurs)
## 1                              ocasional
## 2                           fq_sem_viaje
## 3                          duracion_prom
## 4                              calif_exp
## 5                              volv_casa
## 6                                trabajo
## 7                                escuela
## 8                                cuidado
## 9                                compras
## 10                                  edad
## 11                               niv_edu
## 12                                genero
## 13                                difisi
## 14                                licond
## 15                             estratoCL
## 16                              tot_pers
## 17                              num_moto
## 18                             num_carro
## 19                              num_bici
## 20                          GrandesZones

1.7 Calcul d’une deuxième AFDM après passage des variables peu discriminantes en supplémentaire

res <- FAMD(VariablesTypoUltraMarcheurs, ncp = 6, sup.var = c(3, 4, 17), graph = FALSE) ##sup.var signifie qu'on passe les variables calif_exp (3), duracion_prom (4) et num_moto (17) en supplémentaires. ncp = 6 --> on ne garde que 6 axes, qui expliquent déjà 40% de l'inertie totale.

#Affichage des plans factoriels
plot(res, (choix = c("ind", "var", "quanti", "quali")), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes", palette=palette(c("grey","red","blue")), title = "Factor map of individuals (pedestrians) and categories")

plot(res, choix = c("var"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes", title = "Variable's factor map", palette=palette(c("grey","red","blue"))) 

plot(res, choix = c("quanti"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes") 

plot(res, choix = c("quali"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

#Affichage de la fréquence cumulée des valeurs propres
get_eigenvalue(res)

##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1   3.395012        12.125041                    12.12504
## Dim.2   2.452448         8.758744                    20.88379
## Dim.3   1.846301         6.593932                    27.47772
## Dim.4   1.643386         5.869236                    33.34695
## Dim.5   1.281049         4.575174                    37.92213
## Dim.6   1.220441         4.358717                    42.28084

fviz_screeplot(res)

1.8 Calcul de la CAH (Classification Ascendante Hiérarchique) et graphiques

nbclasses <- 6 # Definition du nombre de classes. 
res.hcpc <- HCPC(res, kk=Inf, nb.clust=nbclasses, consol=FALSE, graph = FALSE) 
#consol = FALSE signifie qu'on n'applique pas a l'issue de la CAH une consolidation par les k-means. graph = FALSE pour ne pas afficher les sorties graphiques. kk=Inf signifie qu'aucune partition prealable n'est réalisée.
plot(res.hcpc, choice="tree", tree.barplot = FALSE, cex = 0.01)

#Affichage dans le plan factoriel des classes
fviz_cluster(res.hcpc,
             shape = 20,
             repel = FALSE,            # Avoid label overlapping
             show.clust.cent = TRUE, # Show cluster centers
             palette = "jco",         # Color palette see ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map",
)

1.9 Création d’un tableau contenant les effectifs par classe et de l’histogramme correspondant

res.hcpc.dataclust <- as.data.frame(res.hcpc[["data.clust"]]) # conversion du tableau obtenu avec la CAH en dataframe
TableauContingence <- data.frame(res.hcpc.dataclust[,"clust"]) # Préparation d'un tableau des effectifs de chaque classe
names(TableauContingence) <- c("clust")
TableauContingence$Count <- 1
TableauContingence <- aggregate(Count ~ clust, TableauContingence, sum)  
View(TableauContingence) # Affichage du tableau d'effectifs de classe
barplot(TableauContingence$Count, names = TableauContingence$clust, main = "Individuals' distribution by cluster") #Affichage des effectifs de chaque classe en histogramme.

1.10 Identification des variables catégorielles ne prenant que deux modalités afin d’alléger la caractérisation des profils.

# Cette opération est nécessaire pour simplifier ensuite les graphiques
TableauVariablesCategorielles <- Filter(is.character, VariablesTypoUltraMarcheurs)

NbModalitesParVariable <- data.frame()

for (i in colnames(TableauVariablesCategorielles)){

  # Recuperation du nombre de modalités pour chaque variable catégorielle
  fichier<- length(unique(TableauVariablesCategorielles[[i]]))
  NbModalitesParVariable <- rbind(NbModalitesParVariable, fichier)
  
}

NbModalitesParVariable$Variable <- colnames(TableauVariablesCategorielles)

VariablesDeuxModalites <- NbModalitesParVariable[which(NbModalitesParVariable[,1] == 2),]
VariablesDeuxModalites <- as.character(VariablesDeuxModalites$Variable)

1.11 Définition du seuil de la valeur V.test à retenir pour conserver les variables décrivant les différentes classes

# Par défaut v.test > |1.96| - Relever le seuil permet de filtrer le nombre de variables qui décriront les classes. On ne garde que celles qui sont les plus significatives
Seuil.v.test <- 3.29
# Avec ce seuil, on accepte un seuil d'erreur à 0,001.
# en savoir plus https://www.medcalc.org/manual/values-of-the-normal-distribution.php

1.12 Création des tableaux décrivant les différentes classes par les variables quanti et quali et affichage des diagrammes

# On ne garde qu'une modalité pour les variables n'ayant que deux modalités
for (i in 1:nbclasses){

  # Récuperation dans des tableaux de la description de chacune des n classes par les variables quantitatives
  b <- as.data.frame(res.hcpc$desc.var$quanti[[i]])
  b <- b[which(colnames(b) == 'v.test')]
  b <- signif(b,3) # pour arrondir les valeurs

  # Récuperation dans des tableaux de la description de chacune des n classes par les variables catégorielles
  c <- as.data.frame(res.hcpc$desc.var$category[[i]])
  d <- as.data.frame(c$label <- rownames(c))
  d <- as.data.frame(c$label <- sub("\\=.*", "", c$label))
  d <- as.data.frame(c[which(c$label %in% VariablesDeuxModalites), ])
  d <- as.data.frame(d[which(d$v.test > 0),])
  e <- as.data.frame(c[which(!c$label %in% VariablesDeuxModalites), ])
  f <- as.data.frame(rbind(d, e))
  g <- as.data.frame(f$label <- NULL)
  c <- signif(f,3) # pour arrondir les valeurs
  c <- c[which(colnames(c) == 'v.test')]

  # Assemblage des tableaux avec les variables quantitatives et qualitatives
  h <- as.data.frame(rbind(b,c))
  h <- as.data.frame(cbind(h, h$label <- row.names(h)))
  colnames(h) <- c("v.test", "label")
  h <- as.data.frame(h[order(h$v.test, decreasing = TRUE), ])
  h <- as.data.frame(h[which(!h$v.test == "Inf"),])
  h <- as.data.frame(h[which(abs(h$v.test) >= Seuil.v.test),])
  
   # Pour tracer l'ensemble des variables significatives sur les diagrammes en barres, il faut au préalable remplacer les V.test infinies (-Inf et Inf). On propose de les remplacer en reprenant les V.test min et max auxquelles on applique un coefficient
  j <- h[which(h$v.test != -Inf),] # création tableau sans les valeurs -Inf
  min.j <- min(j$v.test) # sélection valeur minimale
  h$v.test[which(h$v.test == -Inf)] <- min.j*1.2 # remplacement des -Inf par valeur minimale (multipliée par un coefficient 1.2)
  j <- h[which(h$v.test != Inf),] # création tableau sans les valeurs Inf
  max.j <- max(j$v.test) # sélection valeur maximale
  h$v.test[which(h$v.test == Inf)] <- max.j*1.2 # remplacement des Inf par valeur maximale (multipliée par un coefficient 1.2)
  
  # pour tracer les diagrammes en barres 
  op <- par(oma=c(0,7,0,0))
  barplot((h$v.test), names = row.names(h), border = "white", horiz = TRUE, las = 1, xlim = c(-5-max(abs(h$v.test)), 5 + max(abs(h$v.test))), cex.names=0.8, main = paste("profil",i, sep=""))
  
  
}

1.13 Récupération de la variable clust, et réintégration de la variable id_hog_pers au jeu de données initial

VariablesTypoUltraMarcheurs$profil <- res.hcpc.dataclust$clust
VariablesTypoUltraMarcheurs$id_hog_pers <- row.names(res.hcpc.dataclust)

# utile par la suite pour récupérer les coord. XY du lieu de résidence du ménage
df <- data.frame(x = VariablesTypoUltraMarcheurs$id_hog_pers)
df <- df %>% separate(x, c("id_hog", "id_pers"))
VariablesTypoUltraMarcheurs$id_hog <- df$id_hog

1.14 Exporter le tableur au format .xlsx

write.xlsx(VariablesTypoUltraMarcheurs, "VariablesTypoUltraMarcheursFINAL.xlsx")
#browseURL("VariablesTypoUltraMarcheursFINAL.xlsx")

2 Cartographie des lieux de résidence des marcheurs suivant leur profil

2.1 Cartes des lieux de résidence (représentation sous forme de surfaces lissées)

2.1.1 Importation du tableau contenant les coordonnées des îlots de résidence des marcheurs

CoordXY_Hogar <- read.xlsx("CoordXY_Hogar_Alea.xlsx")

2.1.2 Jointure des coordonnées XY du lieu de résidence dans le tableau des marcheurs

VariablesTypoUltraMarcheurs <- merge(VariablesTypoUltraMarcheurs, CoordXY_Hogar, by = "id_hog")

2.2 Importation de la couche SIG des UTAM (Unités territoriales d’analyse de la mobilité)

UTAM <- st_read("EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3.gpkg")

## Reading layer `EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3' from data source `F:\Rennes\Recherche\Thèses_suivies\Arthur Ducasse\Article_Typo_pietons\EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 132 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N

2.3 Création d’un semis de points (objet sf) à partir des coordonnées XY des lieux de résidence des marcheurs

ResidMarcheurs<- st_as_sf(VariablesTypoUltraMarcheurs, coords=c("X_Hog_Alea","Y_Hog_Alea"), crs = st_crs(UTAM)$srid)
plot(st_geometry(ResidMarcheurs))

2.4 Préparation du jeu de données pour l’élaboration des cartes

# étape nécessaire en vue de créer des cartes lissées donnant à voir la concentration des lieux de résidence des marcheurs suivant leur profil
# pour créer une liste des profils de marcheurs, dans l'ordre croissant des classes
TypoProfilMarcheurs <- VariablesTypoUltraMarcheurs[order(VariablesTypoUltraMarcheurs$profil, decreasing = FALSE), ]
ListProfils <- unique(TypoProfilMarcheurs$profil)
str(TypoProfilMarcheurs)

## 'data.frame':    14218 obs. of  26 variables:
##  $ id_hog       : chr  "10288" "10324" "10496" "10500" ...
##  $ ocasional    : chr  "regular" "regular" "regular" "ocas" ...
##  $ fq_sem_viaje : num  2211 257 780 0 0 ...
##  $ duracion_prom: num  15.8 10 8 30 45 ...
##  $ calif_exp    : num  5 4 4 5 5 3 4 5 5 5 ...
##  $ volv_casa    : num  1005 0 390 130 130 ...
##  $ trabajo      : num  0 0 0 0 0 ...
##  $ escuela      : num  0 257 0 0 0 ...
##  $ cuidado      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ compras      : num  1206 0 0 130 130 ...
##  $ edad         : chr  "sup_65" "men_5" "sup_65" "40-64" ...
##  $ niv_edu      : chr  "univ" "prim_sec" "prim_sec" "univ" ...
##  $ genero       : chr  "muj" "hom" "muj" "hom" ...
##  $ difisi       : chr  "difisi_no" "difisi_no" "difisi_si" "difisi_no" ...
##  $ licond       : chr  "licond_no" "licond_no" "licond_no" "licond_si" ...
##  $ estratoCL    : chr  "s_5-6" "s_5-6" "s_5-6" "s_5-6" ...
##  $ tot_pers     : num  2 4 2 2 2 1 2 2 2 2 ...
##  $ num_moto     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num_carro    : num  0 1 1 4 4 1 1 1 1 1 ...
##  $ num_bici     : num  0 1 0 2 2 0 1 0 1 1 ...
##  $ GrandesZones : chr  "Norte" "Oeste" "Oeste" "Oeste" ...
##  $ profil       : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ id_hog_pers  : chr  "10288-1" "10324-4" "10496-1" "10500-1" ...
##  $ ID_MZA_EMU   : chr  "11001191110709" "11001155060105" "11001163030422" "11001163030422" ...
##  $ X_Hog_Alea   : num  603204 599329 599155 599215 599215 ...
##  $ Y_Hog_Alea   : num  520164 515043 514540 514599 514599 ...

##Création des objets sf dont seront extrait les toponymes qui seront affichés sur les cartes.
#Pour les lieux de résidence
topores1 <- UTAM[which(UTAM$UTAM_NOMBR == "CHAPINERO" | UTAM$UTAM_NOMBR == "USAQUEN"| UTAM$UTAM_NOMBR == "EL PRADO"),]
topores2 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "CARVAJAL" | UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "CHAPINERO"),]
topores3 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "ALAMOS"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO"| UTAM$UTAM_NOMBR == "CHAPINERO"),]
topores4 <- UTAM[which(UTAM$UTAM_NOMBR == "EL LUCERO"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "20 DE JULIO" | UTAM$UTAM_NOMBR == "CIUDAD USME"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO" | UTAM$UTAM_NOMBR == "CHIA" | UTAM$UTAM_NOMBR == "CAJICA" | UTAM$UTAM_NOMBR == "EL ROSAL" | UTAM$UTAM_NOMBR == "BOJACA" | UTAM$UTAM_NOMBR == "FACATATIVA"| UTAM$UTAM_NOMBR == "LA CALERA" | UTAM$UTAM_NOMBR == "SIBATE"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "SUBA"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "FACATATIVA"| UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "SAN ISIDRO - PATIOS"| UTAM$UTAM_NOMBR == "COTA"| UTAM$UTAM_NOMBR == "TABIO"| UTAM$UTAM_NOMBR == "ZIPAQUIRA"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "LA CANDELARIA"),]
topores5 <- UTAM[which(UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "CIUDAD USME" | UTAM$UTAM_NOMBR == "20 DE JULIO"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "LA CALERA"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "EL ROSAL"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "BOJACA"| UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "LA CANDELARIA"),]
topores6 <- UTAM[which(UTAM$UTAM_NOMBR == "EL LUCERO"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "VENECIA"| UTAM$UTAM_NOMBR == "ENGATIVA" | UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "DIANA TURBAY" | UTAM$UTAM_NOMBR == "GRAN YOMASA"),]

#Sauvegarde des objets sf dans une liste
toporesid <- list(topores1, topores2, topores3, topores4, topores5, topores6)

2.5 Production d’une série de cartes montrant la concentration des lieux de résidence des marcheurs selon leur profil

#Définition du nombre de classes
nclass <- 5

#création d'un ensemble de palettes de couleurs (on en crée ici plus qu'il n'en faut, pour avoir de la marge si l'on souhaite augmenter le nombre de classes)
pal1 <- hcl.colors(nclass, "YlOrRd", rev = TRUE)
pal2 <- hcl.colors(nclass, "Greens 3", rev = TRUE)
pal3 <- hcl.colors(nclass, "Reds 3", rev = TRUE)
pal4 <- hcl.colors(nclass, "Grays", rev = TRUE)
pal5 <- hcl.colors(nclass, "Oslo", rev = FALSE)
pal6 <- hcl.colors(nclass, "Red-Purple", rev = TRUE)
pal7 <- hcl.colors(nclass, "Inferno", rev = TRUE)
pal8 <- hcl.colors(nclass, "Reds", rev = TRUE)
pla9 <- hcl.colors(nclass, "Blues 3", rev = TRUE)
pal10 <- hcl.colors(nclass, "Purples 3", rev = TRUE)

#Sauvegarde des palettes dans une liste
pal <- list(pal1, pal2, pal3, pal4, pal5, pal6, pal7, pal8, pla9, pal10)

#Définition des limites des UTAM comme emprise pour le lissage
Emprise <- as.owin(as(st_union(st_buffer(UTAM, 1000)), "Spatial"))

#Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(oma=c(3.5,0,3,0)+0.1, mar = c(0, 0.5, 1.2, 0.5))

plot.new()

#Pour découper la fenêtre en 1 ligne et 6 colonnes (6 profils)
par(mfrow=c(3,2))

#Boucle pour produire et cartographier une surface de densité par profil
for (i in ListProfils){
  
  #Récuperation des jeux de données par profil
  fichier <- TypoProfilMarcheurs[which(TypoProfilMarcheurs$profil == i),]
  
  #Pour récupérer les coordonnées des lieux de résidence
  pts <- st_coordinates(st_geometry(st_as_sf(fichier, coords=c("X_Hog_Alea","Y_Hog_Alea"), crs = st_crs(UTAM)$srid)))
  
  #Pour créer un objet ppp (format spatstat) et y intégrer l'emprise
  fichier.ppp <- ppp(pts[,1], pts[,2], window = Emprise)
  
  # pour définir la taille du rayon (en m.)
  rayon <- 1000
  
  # pour définir la taille du pixel (en m.)
  resol <- 100 # ici 100m x 100m = 1 ha.
  
  # pour calculer la surface de densité (rayon lissage : 1000 m et résolution spatiale de l'image : 100m x 100m = 1 ha)
  cartelissee <- density(fichier.ppp, sigma = rayon, kernel = "gaussian", eps = resol)
  
  #Conversion de la surface lissée au format raster
  cartelissee.raster <- raster(cartelissee)
  crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
  
  #Découpage du raster sur l'emprise des UTAM
  cartelissee.raster <- mask(cartelissee.raster, UTAM) 
  
  #Passage des densités à des effectifs (multiplication des densités par la surface du cercle)
  values(cartelissee.raster) <- values(cartelissee.raster) * pi*rayon**2

  #Définition des seuils de classes (discrétisation en intervalles égaux)
  bks <- seq(from = cellStats(cartelissee.raster, stat = min), 
            to = cellStats(cartelissee.raster, stat = max), 
            by = (cellStats(cartelissee.raster, stat = max) - cellStats(cartelissee.raster, stat = min)) / nclass)
  
  #Reclassification de la surface lissée
  cartelissee.reclass <- cut(cartelissee.raster, breaks = bks, include.lowest = FALSE, right = TRUE, dig.lab = 3, ordered_result = FALSE)
  
  #Vectorisation de la surface reclassée (calcul un peu long)
  cartelissee.vecteur <- as(rasterToPolygons(cartelissee.reclass, n=4, na.rm=TRUE, digits=12, dissolve=TRUE), "sf")
  
  #Tracer la carte
  plot(st_geometry(UTAM), border = "white", bg= "grey90")
  
  typoLayer(
    x = cartelissee.vecteur,
    var="layer",
    col = unlist(pal[as.numeric(i)]),
    lwd = 0.1,
    border = unlist(pal[as.numeric(i)]),
    legend.pos = "n",
    add = TRUE)
  
   labelLayer(
    x = toporesid[[as.numeric(i)]], 
    txt = "UTAM_NOMBR", 
    col= "black", 
    cex = 0.3, 
    font = 4, 
    halo = TRUE, 
    bg = "white", 
    r = 0.05, 
    overlap = FALSE, 
    show.lines = FALSE)
  
  legendChoro(
         pos = "bottomright",
         title.txt = "Number of walker residence\nplaces in a 1 km radius.",
         title.cex = 0.6,
         breaks = bks, 
         values.rnd = 2,
         nodata = FALSE,
         col = unlist(pal[as.numeric(i)]),
         border = "white",
         horiz = FALSE
       )

  title(main =paste("Cluster",i, sep="-"))
  
  mtext(text = paste0("n = ", nrow(fichier), " walkers"), 
      side = 3, line = -2, adj = 0.1, cex = 0.7)
  
}

barscale(
  lwd = 1.5,
  cex = 0.6,
  pos = "bottomleft",
  style = "pretty",
  unit = "m"
)

north(pos = "topright")

#Pour afficher le titre principal et la source
mtext("Concentration of walkers' home location by profile in the metropolitan area of Bogota", cex = 1.3, side=3,line=1,adj=0.5,outer=TRUE)
mtext("   Sources : Union Temporal STEER & CNC, EMD, 2019 - SDM - radius buffer : 1000 m, resolution : 1 ha", side=1, line=1, adj=0, cex=0.6, outer=TRUE)

3 Cartographie des lieux de destination associés au premier déplacement des marcheurs depuis leur domicile

3.1 Importation de la couche SIG représentant les lieux de première destination (centroïdes des ZAT - Zones d’analyse du transport)

#Le premier trajet à pied des marcheurs est répété entre 4 et 5 fois dans la semaine, ce qui le rend assez représentatif des déplacements à pied des marcheurs, puisque le deuxième trajet est un peu moins répété, alors que le 3e est bien plus occasionnel. 
destino <- st_read("ZAT_PrimerDestino.gpkg")

## Reading layer `ZAT_PrimerDestino' from data source 
##   `F:\Rennes\Recherche\Thèses_suivies\Arthur Ducasse\Article_Typo_pietons\ZAT_PrimerDestino.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 14655 features and 31 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 545952 ymin: 453808 xmax: 652779.8 ymax: 600655.7
## Projected CRS: WGS 84 / UTM zone 18N

plot(st_geometry(destino), pch = "+", cex = 0.5)

3.2 Création d’une liste de marcheurs dont la ZAT de destination du premier déplacement est renseignée

listdestino <- as.character(destino$id_concat)

3.3 Sélection des marcheurs ayant une ZAT de destination connue pour leur premier déplacement

TypoProfilMarcheurs <- VariablesTypoUltraMarcheurs[VariablesTypoUltraMarcheurs$id_hog_pers %in% listdestino, ]

3.4 Jointure pour récupérer le profil (clust) sur le lieu de destination

destino <- merge(destino, TypoProfilMarcheurs[c(22:23)], by.x = "id_concat", by.y = "id_hog_pers")

3.5 Construction d’une liste contenant des objets sf dont seront extrait les toponymes à figurer sur les cartes des lieux de destination

##Création des objets sf dont seront extrait les toponymes qui seront affichés sur les cartes

topodest1 <- UTAM[which(UTAM$UTAM_NOMBR == "USAQUEN" | UTAM$UTAM_NOMBR == "CHICO LAGO"),]
topodest2 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "TEUSAQUILLO" | UTAM$UTAM_NOMBR == "CHICO LAGO"| UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "BRITALIA"),]
topodest3 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "TEUSAQUILLO" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "20 DE JULIO" | UTAM$UTAM_NOMBR == "GRAN YOMASA" | UTAM$UTAM_NOMBR == "ENGATIVA" | UTAM$UTAM_NOMBR == "TENJO" | UTAM$UTAM_NOMBR == "VENECIA"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL"| UTAM$UTAM_NOMBR == "SUBA"| UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "FONTIBON"| UTAM$UTAM_NOMBR == "CHICO LAGO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "TUNJUELITO"),]
topodest4 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "FUNZA" | UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO" | UTAM$UTAM_NOMBR == "CHIA" | UTAM$UTAM_NOMBR == "CAJICA" | UTAM$UTAM_NOMBR == "EL ROSAL" | UTAM$UTAM_NOMBR == "BOJACA" | UTAM$UTAM_NOMBR == "FACATATIVA"| UTAM$UTAM_NOMBR == "LA CALERA" | UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "ZIPAQUIRA"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "COTA"| UTAM$UTAM_NOMBR == "SIBATE"| UTAM$UTAM_NOMBR == "TABIO"| UTAM$UTAM_NOMBR == "SUBA"| UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "FONTIBON"| UTAM$UTAM_NOMBR == "VENECIA"| UTAM$UTAM_NOMBR == "20 DE JULIO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "CHAPINERO"),]
topodest5 <- UTAM[which(UTAM$UTAM_NOMBR == "TEUSAQUILLO"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "CIUDAD USME" | UTAM$UTAM_NOMBR == "20 DE JULIO"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "LA CALERA"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "BOJACA"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "SIBATE"| UTAM$UTAM_NOMBR == "ISMAEL PERDOMO"| UTAM$UTAM_NOMBR == "EL ROSAL"| UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "SANTA CECILIA"| UTAM$UTAM_NOMBR == "FONTIBON"| UTAM$UTAM_NOMBR == "TENJO"| UTAM$UTAM_NOMBR == "LOS ALCAZARES"),]
topodest6 <- UTAM[which(UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "SOACHA PUEBLO"| UTAM$UTAM_NOMBR == "VENECIA"| UTAM$UTAM_NOMBR == "TEUSAQUILLO"| UTAM$UTAM_NOMBR == "ENGATIVA"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "KENNEDY CENTRAL"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "RESTREPO"),]

#Sauvegarde des objets sf dans une liste
topodest <- list(topodest1, topodest2, topodest3, topodest4, topodest5, topodest6)

3.6 Production d’une deuxième série de cartes représentant la concentration des lieux associés à la première destination des marcheurs

#Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(oma=c(3.5,0,3,0)+0.1, mar = c(0, 0.5, 1.2, 0.5))

plot.new()

#Pour découper la fenêtre en 2 lignes et 2 colonnes (4 profils)
par(mfrow=c(3,2))

#Définition du nombre de classes
nclass <- 5

#Définition d'une boucle pour produire et cartographier une surface de densité par profil
for (i in ListProfils){
  
  #Récuperation des jeux de données par profil
  fichier <- destino[which(destino$profil == i),]
  
  #Récupération des coordonnées des lieux de résidence
  pts <- st_coordinates(st_geometry(fichier))
  
  #Création d'un objet ppp (format spatstat) et intégration de l'emprise
  fichier.ppp <- ppp(pts[,1], pts[,2], window = Emprise)
  
  #Définition de la taille du rayon (en m.)
  rayon <- 1000
  
  #Définition de la taille du pixel (en m.)
  resol <- 100 # ici 100m x 100m = 1 ha.
  
  #Calcul de la surface de densité (rayon lissage : 1000 m et résolution spatiale de l'image : 100m x 100m = 1 ha)
  cartelissee <- density(fichier.ppp, sigma = rayon, kernel = "gaussian", eps = resol)
  
  #Conversion de la surface lissée au format raster
  cartelissee.raster <- raster(cartelissee)
  crs(cartelissee.raster) <- st_crs(UTAM)$srid
  
  #Découpage du raster sur l'emprise des UTAM
  cartelissee.raster <- mask(cartelissee.raster, UTAM) 
  
  #Passage des densités à des effectifs (multiplication des densités par la surface du cercle)
  values(cartelissee.raster) <- values(cartelissee.raster) * pi*rayon**2

  #Définition des seuils de classes (discrétisation en intervalles égaux)
  bks <- seq(from = cellStats(cartelissee.raster, stat = min), 
            to = cellStats(cartelissee.raster, stat = max), 
            by = (cellStats(cartelissee.raster, stat = max) - cellStats(cartelissee.raster, stat = min)) / nclass)
  
  #Reclassification de la surface lissée
  cartelissee.reclass <- cut(cartelissee.raster, breaks = bks, include.lowest = FALSE, right = TRUE, dig.lab = 3, ordered_result = FALSE)
  
  #Vectorisation de la surface reclassée (calcul un peu long)
  cartelissee.vecteur <- as(rasterToPolygons(cartelissee.reclass, n=4, na.rm=TRUE, digits=12, dissolve=TRUE), "sf")
  
  #Tracé de la carte
  plot(st_geometry(UTAM), border = "white", bg= "grey90")
  
  typoLayer(
   x = cartelissee.vecteur,
   var="layer",
   col = unlist(pal[as.numeric(i)]),
   lwd = 0.1,
   border = unlist(pal[as.numeric(i)]),
   legend.pos = "n",
   add = TRUE)
  
   labelLayer(
    x = topodest[[as.numeric(i)]], 
    txt = "UTAM_NOMBR", 
    col= "black", 
    cex = 0.3, 
    font = 4, 
    halo = TRUE, 
    bg = "white", 
    r = 0.05, 
    overlap = FALSE, 
    show.lines = FALSE)
  
  legendChoro(
         pos = "bottomright",
         title.txt = "Number of destination places\nin a 1 km radius.",
         title.cex = 0.6,
         breaks = bks, 
         values.rnd = 2,
         nodata = FALSE,
         col = unlist(pal[as.numeric(i)]),
         border = "white",
         horiz = FALSE
       )

  title(main =paste("Cluster",i, sep="-"))
  
  mtext(text = paste0("n = ", nrow(fichier), " walkers"), 
      side = 3, line = -2, adj = 0.1, cex = 0.7)
  
}

barscale(
  lwd = 1.5,
  cex = 0.6,
  pos = "bottomleft",
  style = "pretty",
  unit = "m"
)

north(pos = "topright")

# Pour afficher le titre principal et la source
mtext("Concentration of walkers' first destination places in the metropolitan area of Bogota", cex = 1.3, side=3,line=1,adj=0.5,outer=TRUE)
mtext("   Sources : Union Temporal STEER & CNC, EMD, 2019 - SDM - radius buffer : 1000 m, resolution : 1 ha", side=1, line=1, adj=0, cex=0.6, outer=TRUE)